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We present a new frequency-domain phenomenological model of the gravitational-wave signal from 
the inspiral, merger and ringdown of non-precessing (aligned-spin) black-hole binaries. The model 
is calibrated to 19 hybrid effective-one-body-numerical-relativity waveforms up to mass ratios of 
1:18 and black-hole spins of |a/m| ~ 0.85 (0.98 for equal-mass systems). The inspiral part of the 
model consists of an extension of frequency-domain post-Newtonian expressions, using higher-order 
terms fit to the hybrids. The merger-ringdown is based on a phenomenological ansatz that has 
been significantly improved over previous models. The model exhibits mismatches of typically less 
than 1% against all 19 calibration hybrids, and an additional 29 verification hybrids, which provide 
strong evidence that, over the calibration region, the model is sufficiently accurate for all relevant 
gravitational-wave astronomy applications with the Advanced LIGO and Virgo detectors. Beyond 
the calibration region the model produces physically reasonable results, although we recommend 
caution in assuming that any merger-ringdown waveform model is accurate outside its calibration 
region. As an example, we note that an alternative non-precessing model, SEOBNRv2 (calibrated 
up to spins of only 0.5 for unequal-mass systems), exhibits mismatch errors of up to 10% for high 
spins outside its calibration region. We conclude that waveform models would benefit most from a 
larger number of numerical-relativity simulations of high-aligned-spin unequal-mass binaries. 

PACS numbers: 


I. INTRODUCTION 

The first direct gravitational wave (GW) detection is 
anticipated sometime during the operation of the Ad¬ 
vanced Laser Interferometer Gravitational-wave Obser¬ 
vatory (aLIGO) [TH3] and Advanced Virgo (AdV) [4] de¬ 
tectors [5], beginning with aLIGO in 2015. One of the 
most likely sources for the first detections, and a rich 
source of scientific information about both fundamental 
physics and astrophysics [6] are the inspiral and merger 
of binary black hole (BBH) systems. Observations and 
measurements of BBH will rely on accurate theoretical 
models of their GW signal, and the construction of such 
models is currently an active research topic [7]- 

To date most effort has focussed on binaries where the 
spin of each black hole (BH) is either zero, or aligned 
with the binary’s orbital angular momentum. In these 
configurations the orbital plane and spin directions re¬ 
main fixed, and the resulting GW signal is far simpler 
than in generic (precessing) configurations. Recent work 
has suggested that aligned-spin models may allow detec¬ 
tion of most (even precessing) binaries [SUTO]. and also 
that accurate approximate generic models can be con¬ 
structed based on an underlying aligned-spin model m- 

Aligned-spin models that include the two BHs’ inspi¬ 
ral, their merger and the ringdown of the final BH, are 
based on a combination of analytic post-Newtonian (PN) 
and effective-one-body (FOB) methods to describe the 
inspiral, and the calibration of phenomenological merger- 


ringdown models to numerical relativity (NR) simula¬ 
tions. The two classes of models are the phenomeno¬ 
logical (“Phenom”) models 0 IMll, which began as 
phenomenological treatments of both the inspiral and 
merger-ringdown, and FOB models [T5ll28] . which have 
used successively more sophisticated versions of the FOB 
approach to describe the inspiral all the way to merger, 
followed by the smooth connection of a ringdown por¬ 
tion; NR waveforms are used to calibrate unknown FOB 
coefficients and free parameters in the merger-ringdown. 

The original motivation of the Phenom approach was 
to produce an approximate and efficient waveform fam¬ 
ily suitable for GW searches (the models are written 
as closed-form analytic expressions in the frequency do¬ 
main), and indeed this practical approach allowed the 
construction of the first aligned-spin model, often re¬ 
ferred to as “PhenomB” [8]. Although some aspects of 
the model were made more accurate in the succeeding 
“PhenomC” model mi, the Phenom approach is still re¬ 
garded by many as approximate, and in particular not 
suitable for parameter estimation. This perception has 
been reinforced by the limited region of parameter space 
over which the aligned-spin PhenomC model was cali¬ 
brated — up to binary mass ratios of only 1:4 (spinning 
up to 1:3), and BH spins of only ajm ^ 0.75 (0.85 for 
equal-mass systems). In this work (and its companion 
article, which we will refer to as Paper 1), we show that 
the phenomenological approach is capable of describing 
BBH waveforms with a high degree of physical fidelity. 
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well within the requirements of aLIGO and AdV, and we 
construct a model that is calibrated to the largest region 
of parameter space to date — up to mass ratios of 1:18, 
and spins up to ajm ^ 0.85 (0.98 for equal-mass sys¬ 
tems). This constitutes the main purpose of this paper, 
to present our new “PhenomD” model, and demonstrate 
its accuracy. 

In contrast, the most recent EOB-NR (SEOB- 
NRv2) [26] model is calibrated to NR waveforms up to 
mass ratio 1:8, and spins up to a/m ^ 0.5. It has been 
shown to be extremely accurate within its calibration re¬ 
gion, and it also appears to produce physically reason¬ 
able waveforms over the full range of BH spins, and up 
to much higher mass ratios [29|.. In this work, how¬ 
ever, we find that the SEOBNRv2 model may not ac¬ 
curately describe the merger-ringdown regime for high 
spins ajm > 0.7. This finding motivates a second pur¬ 
pose of this paper: to make clear that the accuracy of 
any merger-ringdown model, Phenom, EOB-NR, or oth¬ 
erwise, is only as good as its NR calibration region. The 
model may give physically plausible results, but its ac¬ 
curacy cannot be guaranteed until it has been checked 
against fully general relativistic NR calculations, and its 
accuracy may well be poor until it has been calibrated 
to those simulations. This seemingly obvious observation 
bears emphasising. It also motivates efforts to quantify 
the accuracy of PN and EOB calculations increasingly 
far back into the inspiral [21 [30]. 

Another important contribution of the Phenom pro¬ 
gramme has been to isolate which combinations of phys¬ 
ical binary parameters will be measurable in GW obser¬ 
vations. Eor example, the previous aligned-spin Phenom 
models Emu exploited the observation that the domi¬ 
nant spin effect on the GW phase is due to a weighted 
combination of the individual BH spins, and the models 
depend on only two physical parameters, the symmetric 
mass ratio and this single effective spin parameter. The 
identification of a simple combination of the in-plane spin 
components in generic binaries m in turn led to a sim¬ 
ple extension of PhenomG to produce a generic-binary 
model, PhenomP [5^ . 

A corollary of this parameter-space reduction is that 
individual spins are expected to be difficult to mea¬ 
sure from GW observations, even if we have a two-spin 
model to hand. Based on previous studies pilHi). and 
a forthcoming study that illustrates in detail the diffi¬ 
culty of measuring individual spins with an aligned-spin 
model m, we also use an effective reduced spin pa¬ 
rameter in certain parts of the PhenomD model. We 
will nonetheless pursue the extension of the Phenom ap¬ 
proach to two spins in future work. 

An additional feature of the PhenomD model is its 
modularity. The separate inspiral and merger-ringdown 
parts of the model are connected by the requirement of 
continuity in the phase and amplitude. This simple con¬ 
struction makes it straightforward to improve and change 
either part of the model independently. We make use of 
this feature to compare versions with alternative choices 


for the inspiral part of the model. 

This paper is organized as follows. In Paper I we dis¬ 
cussed in detail the numerical simulations we have used, 
and in particular presented studies of the accuracy of the 
new NR waveforms that we have produced. In this paper 
we re-visit these waveforms, but from the point of view 
of GW applications, and assess their accuracy in terms of 
their noise-weighted inner product (match). The match 
is defined in Sec. m along with techniques that we use 
to estimate the match between NR waveforms over fre¬ 
quency ranges that extend beyond those where we have 
NR data. In Sec. imi we summarize the waveforms that 
we use, and present our match-based accuracy analysis. 
In Secs. [V| and [VT| we give details of procedure we use to 
construct our models of the signal phase and amplitude, 
over three frequency regions. More details are provided 
in Paper I, but here we summarize the approach, and 
its use across all of the waveforms used to calibrate our 
model, and the accuracy of the final models for each of 
its six constituent parts (th ree phase parts, and three 
amplitude parts). In Sec.j^we assess the final complete 
model’s accuracy by calculating matches against both the 
waveforms used for calibration, and an additional set of 
waveforms that were not used for calibration. We discuss 
the accuracy of our single-reduced-spin approximation, 
and our choice of the minimal set of waveforms neces¬ 
sary for the model. In Sec. we compare against the 
SEOBNRv2 model, illustrating the high-spin, unequal- 
mass region where we find disagreement between the two 
models; this is outside the calibration region of SEOB- 
NRv2. In Appendix we revisit the agreement between 
our new model and the original NR data by transforming 
PhenomD to the time domain, and in Appendixj^we list 
the PN inspiral coefficients used in our model. 


II. PRELIMINARIES 

A. Outline of the model 


We describe a BBH system by the following param¬ 
eters. The masses are rui and m 2 , where we choose 
mi > m 2 , and the total mass is M = mi + m 2 . The 
mass ratio of the binary is denoted g = mi/m 2 > I, and 
the symmetric mass ratio is 7 ^ = mim 2 /M^. The BH spin 
angular momenta are Si and S 2 which we assume to be 
parallel to the direction of the orbital angular momen¬ 
tum, L. In this work we restrict ourselves to aligned-spin 
(non-precessing) systems, and so are only concerned with 
the dimensionless spin parameters defined as 

Si’L 


with Xi e [-1,1]- 

In previous aligned-spin Phenom models, we have pa¬ 
rameterized the spin depedence of the model by a single 
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effective spin parameter uni], 


Xeff = 


mixi +m2X2 

M 


( 2 ) 


This is based on the observation that it is a weighted 
sum of the spins that constitutes the dominant effect of 
the spin on the inspiral of the binary. In PN theory, 
the leading-order spin effect on the binary’s phasing is in 
fact [Ml [^[37], 


3877 . . 

Xpn = Xeff - Yi3 ^ 

and we have seen evidence in previous work that this 
is in general a better parameter to use also in inspiral- 
merger-ringdown (IMR) models [Ml- In this work we use 
a combination of spin parametrizations. The phenomeno¬ 
logical calibrations to NR waveforms are parameterized 
by Xpn, normalized such that its range is from -1 to 1 
for all mass ratios, 


^ 1-7677/113' 

The final BH is correctly parameterized by the final 
mass Mf and spin a/, and for this reason the final mass 
and spin estimates that we use (see Paper 1), are parame¬ 
terized by a different spin combination. Si ^ 82 - Finally, 
our inspiral model is based on the standard frequency- 
domain PN approximant, “TaylorF2” [MHlOl, and this is 
parameterized by both spins, xi and X 2 - The final result 
is a model that depends on both spins xi and X 2 , but the 
calibration to hybrid EOB+NR waveforms is parameter¬ 
ized by different combinations of xi and X 2 for the in¬ 
spiral, merger and ringdown parts of the model. Most of 
the hybrid waveforms are for equal-spin x = Xi = X 2 sys¬ 
tems, so we can guarantee our model’s accuracy only for 
these configurations. However, as we discuss in Sec. |IXB , 
the X approximation is extremely accurate for most re¬ 
gions of parameter space, and in those where it is not 
(higher mass ratios and high parallel spins), the innacur- 
racy is unlikely to have any influence on GW astronomy 
applications with aLIGO or AdV. 

The PhenomD model provides expressions for the £ = 
2, |m| = 2 spin-weighted spherical-harmonic modes of 
the GW signal, since these are the dominant modes in 
aligned-spin systems. The full signal as a function of the 
physical parameters E G (M, 77 , Xi, X 2 ) and the observer’s 
orientation (^, 0) with respect to the orbital angular mo¬ 
mentum of the binary, is given by, 

^(/; 2 , 0,4>) = h+{f] 5,6», 4)) - ihx if; s, 0, f) (5) 

= E h2mif;E) -^Y2mi0,4>), (6) 

m=-2,2 

where /i 2 ,- 2 (/) = ^2 2 (“/)- We express h 22 {f) in terms 
of the signal amplitude and phase by 

/i22(/;S)=^(/;S)e-W;H)^ (7) 


and it is models of A{f; 5) and 0(/; 5) that we provide. 
Note also that the total mass M provides an overall scale 
for our waveforms, so the physical parameters over which 
the model has been explicitly constructed are 77 , xi and 
X 2 (with the spins treated in combinations as described 
above). 

As ingredients in our model construction, we use hy¬ 
brid waveforms, where the early inspiral is described 
by the nn-c alibrated SEOBv2 model (see Paper 1, 
and Sec. IV below), and the late inspiral and merger- 
ringdown by NR waveforms. The mass and spin of the 
final BHs, Mf and aj, which are key parts of the merger- 
ringdown model, are provided by fits to NR data. The 
details of the hybrid construction, and of the final mass 
and final spin fits, are given in Paper 1. 

We model separately three frequency regimes of the 
waveform. The first region covers the inspiral, up to the 
frequency Mf = 0.018. Here the information is pre¬ 
dominantly from the analytical EOB inspiral waveforms, 
although there is some information at higher frequencies 
from the early parts of the longer NR waveforms; the fre¬ 
quency at which each hybrid switches to an NR waveform 
is provided in Tab.|Tj The second two regions are informed 
purely from NR data. We note that in principle one 
could also construct the individual inspiral and merger- 
ringdown models separately from PN or EOB models (for 
the inspiral) and NR data (for the merger-ringdown), 
without constructing any hybrid waveforms. In this work 
we chose to use hybrid waveforms, because they allow us 
to use the maximum NR information (which influences 
to some extent our inspiral model), and allows for a con¬ 
sistent choice of calibration points in parameter space for 
both the inspiral and merger-ringdown. 

The resulting model is modular: we are free to use a 
different inspiral model, or a different merger-ringdown 
model, as we wish. This introduces a flexibility that 
was not present in previous models. If in the future we 
have access to a more accurate inspiral model (EOB, 
PN, or otherwise), or more accurate merger-ringdown 
model (e.g., calibrated to waveforms over a larger re¬ 
gion of parameter space), then we can easily replace that 
part of the model without any additional tuning. The 
model calculates appropriate time- and phase-shifts (a 
linear correction to the frequency-domain phase) to en¬ 
sure that the phase connects smoothly between the in¬ 
spiral and merger-ringdown, and the model of the am¬ 
plitude in the intermediate region between inspiral and 
merger-ringdown is constructed such that the function is 
continuous. 


B. Matches 

To assess the accuracy of our model and generally 
quantify the (dis)agreement between two waveforms hi 
and h 2 (real-valued in the time domain), we use the stan¬ 
dard inner product weighted by the power spectral den- 
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sity of the detector 5'n(/)- It is defined as [55] . 


{hi, ^ 2 ) = 4 Re 



hi{f) Kif) 

Snif) 


( 8 ) 


The match between two waveforms is defined as the 
inner product between normalised waveforms {h = 
h /maximised over relative time and phase shifts 
between the two waveforms, 


M{hi, h 2 ) = max (hi, ^ 2 )- (9) 


A time- and phase-shift has no significance for the phys¬ 
ical fidelity of an aligned-spin waveform — they corre¬ 
spond, respectively, to a change in the merger time of 
the binary, and of the initial phase of the binary, i.e., an 
overall rotation. 

Results will be quoted in terms of the mismatch A4, 
defined as. 


M(hi,h2) = l-M(hi,h2). ( 10 ) 


We use two noise spectra in this work: the “early 
aLIGO” spectrum, which approximates the detector re¬ 
sponse during the first observing run, planned for late 
2015, and the “zero-detuned high-power” (zdethp) spec¬ 
trum, which is the design goal of aLIGO that is antici¬ 
pated by 2019-20 [41]. Galculations with the early aLIGO 
curve use a lower cutoff frequency of /min = 30 Hz, and 
zdethp calculations are carried out with /min = 10 Hz. 
In both cases, we use /max = 8000 Hz which is greater 
than the highest frequencies contained in the signals we 
are considering. 

In various steps of the model construction in this pa¬ 
per, we are interested in analyzing the agreement of wave¬ 
form sections that are only defined over a certain fre¬ 
quency range. (A good example are NR waveforms that 
are typically too short to fill the entire aLIGO frequency 
band.) In these cases, one could reduce the integration 
limits in ^ to the frequency range defined by the wave¬ 
form sections, but the resulting matches would be dif¬ 
ficult to interpret as they have no direct application in 
GW searches. Here instead, we ask the question “What 
influence does the difference in a certain signal part have 
on the full waveform, assuming all other parts are per¬ 
fectly modeled?” We address this question by aligning 
the signal parts that we wish to compare as if they were 
hybridized with a common model of the remaining signal 
and set the phase difference for this particular alignment 
to zero over all frequencies that are not covered by the 
waveform sections we consider. To construct the full in¬ 
tegrand in (§, we additionally need a model of the am¬ 
plitude, which we take from our final PhenomD model, 
although this particular choice is far less important than 
the phase disagreement we wish to quantify. We can then 
use a standard algorithm to calculate the mismatch be¬ 
tween both signals, and due to their simple form in the 
frequency domain, time and phase shifts will be properly 
taken into account across the entire signal. More details 
and a full discussion of this approach is given in m- 


III. NUMERICAL-RELATIVITY WAVEFORMS 


We calibrated the PhenomD model with publicly avail¬ 
able NR waveforms from the Simulating Extreme Space- 
times (SXS) collaboration [43|, and a set of new simula¬ 
tions produced with the BAM code O ES]. Details of 
the new BAM simulations and their numerical accuracy 
are presented in Paper I. Here we summarize the 19 NR 
waveforms that we used to calibrate the model, and the 
additional waveforms that were used to further test its 
accuracy. 

Our two main goals are to extend the parameter- 
space coverage of aligned-spin phenomenological models 
to higher mass ratios, and to improve the overall accuracy 
to well within the requirements of GW detection and pa¬ 
rameter estimation with Advanced LIGO and Virgo; in 
practice we consider a mismatch error of less than 1 % to 
be sufficient. The first goal dictated our choice of new 
NR simulations. 

The previous aligned-spin phenomenological models, 
PhenomB [ 8 ] and PhenomC mi, were constructed from 
waveforms up to mass ratios of 1:4, and (equal) spins up 
to ±0.75 (with ±0.85 for equal-mass binaries), although 
spinning-binary waveforms were used only up to mass- 
ratio 1:3. We found in constructing those models that it 
was sufficient to use only four or five NR waveforms in 
each direction of parameter space. This suggests that we 
can construct a model across the entire ( 77 , y) parameter 
space with only 25 waveforms. 

Five waveforms equally spaced in 77 would be placed 
at 77 = (0.25, 0.20,0.15, O.IO, 0.05). (In the current model 
we do not include extreme-mass-ratio 77 ^ 0 waveforms, 
e.g.. Refs. [46l [T7] . but we plan to use these to complete 
our parameter-space mapping in future work). We focus 
on simulations at mass ratios q = 1,4,8,18, which cor¬ 
respond to 77 ^ (0.25, 0.16, 0.10, 0.05); we find that wave¬ 
forms at 77 ~ 0.2 are not necessary to produce an accurate 
model, although the model is tested against waveforms 
at g = 2, 3 (77 = 0.222, 0.1875). 

We produced new waveforms with the BAM code up to 
mass ratio 1:18, and for a range of spins. At lower mass 
ratios we have also used publicly available waveforms, 
which were produced by the SXS collaboration using the 
Spectral Einstein Gode (SpEG). In particular, their cat¬ 
alogue provides waveforms for equal-mass binaries with 
high BH spins of —0.95 and ±0.98. The parameter space 
coverage of NR waveforms used in previous models, and 
in our new model, are shown in Fig. and the details 
of the waveforms that we used are summarized in Tab. [H 
We tested the model against an extended sent of wave¬ 
forms, and this is described in more detail in Sec. |IX| and 

Tab.lini 

The accuracy of the new BAM simulations was dis¬ 
cussed in some detail in Paper I. In this work we are 
interested in constructing accurate waveform models for 
GW astronomy with aLIGO and AdV. In that context, 
an important accuracy measure is the mismatch between 
the waveforms with respect to the aLIGO noise spectrum. 
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FIG. 1: Parameter space over which the PhenomD model 
has been calibrated. The locations in parameter space of the 
calibration waveforms are indicated by red points. Also shown 
are the calibration points for the SEOBNRv2 (green) and 
PhenomC (blue) models. 


We calculate the mismatch between the numerical wave¬ 
forms following the procedure outlined in Sec. El in 
particular, we take into account the inspiral signal power, 
allowing us to calculate mismatches for low-mass sys¬ 
tems, and reliably infer the (typically larger) mismatches 
in these systems due to any errors in the merger-ringdown 
waveforms. This procedure tends to estimate larger mis¬ 
matches than integrating Eq. (§ over only the frequency 
range of the NR waveforms, as in, e.g.. Ref. [48], and is a 
more conservative estimate of the mismatch error in the 
NR waveforms. 

We consider the effect of two sources of error on the 
mismatch: the errors due to (1) finite numerical reso¬ 
lution, and (2) finite waveform extraction radius. In all 
cases we have found the overall mismatch error from these 
sources to be < 0.5%. Here we focus on two configura¬ 
tions, g = 4, xi = X 2 = X = 0.75 (AlO), and nonspinning 
^ = 18 (A18). 

Fig-i shows the mismatch error due to numerical res¬ 
olution. In the q = A configuration, the reference sim¬ 
ulation uses a base grid size of 112^ points, with the 
finest grid spacing being h^in = M/230. Comparisons 
are made against simulations with the same resolution, 
but a base grid size of 96^ points, and an 80^ simulation 
with the resolution scaled to give the same physical grid 
sizes as in the 96^ simulation. The solid black line shows 
the mismatch between the 112-point and 96-point simu¬ 
lations, i.e., simulations where only the physical grid sizes 
were changed. This change introduces a mismatch error 
of at most ^ 0.01%. The dashed black line shows the mis¬ 
match between the 112-point and 80-point simulations, 
i.e., both the physical grid sizes and the numerical reso¬ 
lution have been reduced. Here the mismatch difference 
is at most ^ 0.1%. 

The orange lines show the mismatch between the 
g = 18 waveforms, with grid sizes of 96^, 120^ and 144^ 
points. These three simulations constitute a convergence 
series, and we have shown in Paper 1 that they exhibit 



PIG. 2: Mismatch error due to numerical resolution, for the 
q — — X — 0.75 (black lines) and non-spinnning 

g = 18 simulations (orange lines). The solid black line shows 
the mismatch between waveform g = 4 112- and 96-point 
simuations, and the dashed black line shows the mismatch 
between the 96- and 80-point simulations. Eor the g = 18 
configuration, the solid orange line shows the mismatch be¬ 
tween the 144- and 120-point simulations, and the dashed or¬ 
ange line shows the mismatch between the 144- and 96-point 
simulations (see text). 



EIG. 3: Mismatch errors due to finite-radius waveform ex¬ 
traction for the 120-point simulations of the same g = 4 case 
as in Eig. Mismatches are between the Rex = 100 M 
waveform and those extracted at Rex = {50, 60, 70,80,90} M 
(from top to bottom). 


evidence of sixth-order convergence. The solid orange 
line shows the mismatch between the 144^ and 120^ sim¬ 
ulations, and the dashed orange line shows the mismatch 
between the 144^ and 96^ simulations. The higher mis¬ 
matches at high mass, compared to the g = 4 configura¬ 
tion, suggests that the merger-ringdown errors are larger 
in this case, although their effect on the mismatches at 
lower masses is comparable. We again conclude that the 
waveforms are accurate to well within our 1% criterion. 

Although the convergence of our simulations is in gen¬ 
eral unclear, we typically find that our 80-point simula¬ 
tions are not in the convergence regime, and are much less 
accurate than higher-resolution simulations. We there¬ 
fore expect that if the mismatch between the 112-point 
and 80-point simulations is no larger than 0.1%, then the 
mismatch between the 96- or 112-point simulations and 
the continuum limit will be lower than this; it will cer¬ 
tainly be lower than the 1% accuracy requirement that 
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we place on our model. 

Fig. i shows the mismatch between waveforms ex¬ 
tracted at different radii. The waveforms were extracted 
at Rex = {50, 60, 70, 80, 90,100} M, and the mismatch 
calculations are performed against the Rex = 100 M 
waveforms. We expect the error to fall of as ^ IjRex^ and 
in general we observe this for our simulations, but only 
for Rex ^ 60 M. Since even the Rex = 50 M waveform 
has a mismatch of only ^ 0.3% with the Rex = 100 M 
waveform, and assuming a 1 /Rex fall-off in waveform ex¬ 
traction error, we expect that the contribution of this 
error to the Rex = 100 M waveforms is less than 0.1%. 

Based on this analysis, we conclude that our simula¬ 
tions are well within the accuracy requirements to con¬ 
struct a waveform model with an overall mismatch error 
of < 1%. 


IV. CHOICE OF INSPIRAL APPROXIMANT 

The early, gradual inspiral of compact binaries and the 
GWs they emit can be accurately modeled by expanding 
the energy and flux of the system into a PN series. De¬ 
pending on how the underlying equations are formulated 
and solved, there is a variety of PN approximants, each 
consistent with the others when truncated at the same 
expansion order. However, as every approximant is for¬ 
mulated with different, mostly implicit, assumptions of 
how higher order terms are treated, the GW signals they 
predict can differ considerably, especially towards higher 
mass ratios, increased spin magnitudes and for increasing 
orbital frequencies mimsi]. There are sophisticated 
methods that aim to improve the convergence and ac¬ 
curacy of PN-based approximants, and one of the most 
successful approaches is the mapping to an EOB system 
[55H57]. 

In the construction of a complete waveform model we 
face the following two issues. First, we need to pick one 
approximant that, to our current knowledge, models the 
inspiral most accrately. Second, this inspiral description 
has to be complemented by NR-based information about 
the merger and ringdown. We briefly summarize our 
strategy to address both issues below and give references 
to the following sections that describe our reasoning in 
more detail. 

Recent studies have indicated that among the family of 
non-precessing inspiral approximants, the EOB approxi¬ 
mant by Tarachini et al. [26] shows the most consistent 
agreement with NR simulations within the calibration 
range of the model [2aEni. In Paper 1, we have per¬ 
formed an independent consistency test between inspiral 
approximants and our set of NR data and confirmed this 
conclusion. (Note that the most recently calibrated ver¬ 
sion of a non-precessing EOB model [28] has not yet been 
included in any of these tests.) Hence, we used the Tara¬ 
chini et al. model (dubbed SEOBNRv2 in the publicly 
available LIGO software library [58]) as our target inspi¬ 
ral approximant, albeit in its original, uncalibrated form 


that does not include NR fitted corrections (we refer to 
this form as as SEOBv2). Specifically, this involves cal¬ 
culating the SEOBNRv2 waveforms with all of the NR 
calibration terms set to zero, to provide an “uncalibrated” 
SEOBv2 calculation of the inspiral waveform. 

We do so because our goal is to explore an alterna¬ 
tive modeling approach that is independent of previous 
NR-informed EOB tuning. In particular, we performed 
dedicated NR simulations outside the calibration range of 
SEOBNRv2, and instead of inheriting higher-order cor¬ 
rections that were fitted in a smaller parameter space re¬ 
gion, we prefer to use the uncalibrated EOB model purely 
in the inspiral regime and hybridize it with NR data of 
the merger and ringdown. 

We are naturally limited by the lengths of the NR 
waveforms, which are different for every simulation. Pre¬ 
vious studies of NR waveform length requirements have 
suggested that PN inspiral waveforms up to 5-10 orbits 
before merger are sufficiently accurate for detection pur¬ 
poses [121 [51]; many more orbits are needed to fulfil more 
stringent accuracy requirements [SH [531 EOl EB , espe¬ 
cially in the high-mass-ratio and high-spin regime that 
we are covering. Many of our NR waveforms are too 
short to allow that. However, previous studies estimated 
the accuracy of PN approximants based on the differences 
between all available approximants at 3.5PN order (with 
highest spin corrections at 2.5PN order at that time). 
One might argue that the EOB approach is more accu¬ 
rate, and therefore comparisons between PN waveforms 
exaggerate the uncertainty in our best current models. 
On the other hand, without fully general-relativistic re¬ 
sults to compare to, one might be sceptical of good agree¬ 
ments between alternative EOB waveforms that are very 
similar by construction. 

Nevertheless, given that we can join EOB with our NR 
data in a much more robust manner than any of the PN 
approximants (see Sec. H of Paper 1 for our full anal¬ 
ysis), we trust that they provide a reasonably accurate 
description of the inspiral up to the point where NR data 
take over. At what frequency this switch from EOB to 
NR happens depends on the length of the individual NR 
simulations. We note that the lowest common starting 
frequency of our NR waveforms is Mf ^0.018, and this 
is where we begin our phenomenological merger-ringdown 
model. Note, however, that our hybridization procedure 
ensures that the maximum amount of NR information is 
used in every point of the parameter space to inform both 
the inspiral and merger-ringdown part of our model. 


V. MODEL OF THE NR REGIME (REGION II) 

We model separately three frequency regions of the 
waveforms. These are indicated in Fig. Region I is 
defined to be the portion of the hybrid that contains the 
optimal blend of NR and SEOBv2 data. Region II is 
the portion of the hybrid that contains purely NR data 
and corresponds to frequencies Mf > 0.018. This re- 
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# 

Code/ID 

q 

g 

Xi 

X2 

X 

Mf 


M/rd M/hyb 

A^gw,nr 

AI 

SXS:BBH:0I56 

I. 

0.25 

-0.95 

-0.95 

-0.95 

0.9681 

0.3757 

0.0713 

0.00522 

22 

A 2 

SXS:BBH:0I5I 

I. 

0.25 

- 0.6 

- 0.6 

- 0.6 

0.9638 

0.4942 

0.0764 

0.00517 

26 

A3 

SXSiBBHiOOOI 

I. 

0.25 

0 . 

0 . 

0 . 

0.9516 

0.6865 

0.0881 

0.00398 

54 

A4 

SXS:BBH:0I52 

I. 

0.25 

0.6 

0.6 

0.6 

0.9269 

0.8578 

0.1083 

0.00501 

42 

A5 

SXS:BBH:0I72 

I. 

0.25 

0.98 

0.98 

0.98 

0.8892 

0.9470 

0.1328 

0.00497 

48 

A 6 

BAM 

4. 

0.16 

-0.75 

-0.75 

-0.75 

0.9846 

0.0494 

0.0614 

0.00713 

15 

A7 

BAM 

4. 

0.16 

-0.5 

-0.5 

-0.5 

0.9831 

0.1935 

0.0649 

0.00716 

18 

A 8 

SXS:BBH:0I67 

4. 

0.16 

0 . 

0 . 

0 . 

0.9779 

0.4715 

0.0743 

0.00665 

28 

A9 

BAM 

4. 

0.16 

0.5 

0.5 

0.5 

0.9674 

0.7377 

0.0906 

0.008II 

26 

AlO 

BAM 

4. 

0.16 

0.75 

0.75 

0.75 

0.9573 

0.8628 

0.1054 

0.00818 

30 

All 

BAM 

8 . 

0.099 

-0.85 

-0.85 

-0.85 

0.9898 

-0.3200 

0.0546 

0.00918 

8 

A 12 

SXS:BBH:0064 

8 . 

0.099 

-0.5 

0 . 

-0.458 

0.9923 

-0.0526 

0.0589 

0.00632 

36 

A13 

SXS:BBH:0063 

8 . 

0.099 

0 . 

0 . 

0 . 

0.9894 

0.3067 

0.0677 

0.00623 

49 

A14 

SXS:BBH:0065 

8 . 

0.099 

0.5 

0 . 

0.458 

0.9846 

0.6574 

0.0838 

0.00615 

66 

A15 

BAM 

8 . 

0.099 

0.85 

0.85 

0.85 

0.9746 

0.8948 

0.1087 

0.01580 

15 

A16 

BAM 

18. 

0.05 

- 0.8 

0 . 

-0.77 

0.9966 

-0.53II 

0.0514 

0.01035 

14 

A17 

BAM 

18. 

0.05 

-0.4 

0 . 

-0.385 

0.9966 

-0.1877 

0.0563 

0.01283 

15 

A18 

BAM 

18. 

0.05 

0 . 

0 . 

0 . 

0.9959 

0.1633 

0.0633 

0.01284 

13 

A19 

BAM 

18. 

0.05 

0.4 

0 . 

0.385 

0.9943 

0.5046 

0.0745 

0.00916 

23 


TABLE I: Hybrid waveform configurations used to calibrate the PhenomD model. For each configuration we list both the mass 
ratio q and symmetric mass ratio 77 , along with the spins xi cind X 2 and the reduced-spin combination, x, which follows from 
Eq. Q. The final BH has mass Mf and dimensionless spin af, and the ringdown signal has frequency M/rd- The frequency 
M/hyb marks the midpoint of the transition region between SEOBv2 inspiral and NR data. The approximate number of NR 
GW cycles in each hybrid is given by iVGw,NR. 


gion is further sub-divided into two regions, Regions Ila 
and Ilb. These divisions correspond to the intermediate 
and merger-ring down models for both the amplitude and 
phase. 

The figures indicate both the frequency ranges over 
which the three parts are connected, but also the ranges 
that are used to calibrate the model’s coefficients to the 
hybrid data. These regions are in general slightly larger 
than those used when piecing together the final model. 

We will refer to other features of these figures in the 
forthcoming sections. 


A. Prom PhenomC to PhenomD 


The merger-ringdown portion of the phase was mod¬ 
elled in PhenomC M using the ansatz. 


V’l>ii(/) = ^ {aif + 0.2! ^ 

+<^3/ + 0^4 + • 


( 11 ) 


The phase was fit over the frequency range [0.1, 1]/rd- 
The reference phase and time of the fit are given by the 
coefficients 0^4 and a^. At the ringdown frequency /rd 
the phase was smoothly connected to a linear function, 
'^rd(/) = a + /^ 2 /, using a tanh transition function. 


We now aim to model the merger-ringdown phase of 
the NR waveforms only from Mf = 0.018, to ensure that 
we include only NR information in this part of the model. 
Fig.j^shows the derivative of the frequency-domain phase 
for the configuration q = xi = X 2 = —0.95. The 
dashed line shows a fit to the phase using the proce¬ 
dure described above; beyond the ringdown frequency 
^/rd = 0.071 the derivative of the phase is constant, 
and in this example the transition is only piecewise con¬ 
tinuous. We see that, while Eq. is able to accurately 
reproduce the phase up to the ringdown frequency, the 
linear approximation at higher frequencies is crude. 

The solid line in Fig. shows a fit to the phase follow¬ 
ing the procedure we use to construct PhenomD, which 
was motivated in detail in Paper 1, and is also described 
in Sec. |VBT] below. This accurately reproduces the main 
features of the phase derivative in the vicinity of the ring- 
down frequency. There is some disagreement at higher 
frequencies, but we note that the accuracy of the NR data 
typically degrades at these frequencies, and the true be¬ 
haviour of (/)'{/) is not clear. 

In the next section we describe the methodology used 
to produce models of the phase and amplitude for the late 
inspiral, merger and ringdown parts of the waveform, i.e., 
those frquencies for which we have NR data. These we 
have denoted Region II; see Fig. We assume that we 
have a valid inspiral approximant that we can join to 
our NR-based Region II model to construct a full IMR 











Mf 

FIG. 4: Phase derivative —0'(/) = —d(l){f)/df (upper panel) 
and amplitude (lower panel) for the q = 1, xi — X‘2 — —0.95 
configuration. The frequency ranges that were used in the 
fits for each section are shown as black double-ended arrows. 
For reference, the frequency Mf = 0.018 is marked with a 
black dashed line. Shaded regions illustrate the boundaries 
between the different regions when constructing the full IMR 
waveform. The ringdown frequency for this case is Mf = 
0.071. 



FIG. 5: Phase derivative 0^(/) for the q = 1, Xi — X 2 = 
—0.95 configuration. The numerical data (dotted) show a 
distinctive extremum at the ringdown frequency, M/rd = 
0.071, indicated by a vertical dashed line. A fit that follows 
an approach similar to that used for PhenomG (dashed) is 
only a crude approximation to the phase for / > /rd , whereas 
the approach used for the PhenomD model (solid) accurately 
models the phase into the ringdown. 


phase, d(\)/df = 0^(/). For this reason we first model 0', 
and then integrate the resulting expression to produce 
the final phase model. We also note that the overall I /77 
dependence in the inspiral, Eq. (27), also holds for the 
merger and ringdown, and so all of our primary fits are 
to r](\)'. 


1. Region Ilb - merger-ring down 


waveform model. The construction of a suitable inspiral 


model (Region I) is given in Sec. VI 


Our current construction requires that the starting fre¬ 
quency of the Region II model must be consistent for all 
waveforms. This imposes the constraint that the starting 
frequency of the NR-based Region II model is the low¬ 
est common GW frequency for which we have NR data, 
Mf r\j 0.018. This is purely based on the available NR 
data and could in principle be pushed towards lower fre¬ 
quencies given longer waveforms. 


B. Phase 

To produce a robust model there are two key require¬ 
ments: ( 1 ) the ansatz must fit the data well, i.e., the fits 
have small residuals to the data, and ( 2 ) the choice of 
ansatz should ideally be chosen to in such a way that the 
coefficients vary smoothly across the parameter space, to 
enable an accurate parameter-space fit in the final model. 

We find that a simple approach is to split Region II 
into an intermediate (Region Ila) and merger-ringdown 
(Region lib) part, and model them separately, as shown 
in Fig. 

The detailed features of the phase through Region II 
are most apparent when we consider the derivative of the 


An example of the derivative of the phase, 0' is shown 
in Fig. I^for a binary with ^ = 1, Xi = X 2 = —0.95. As 
described in Paper 1, we propose the following ansatz to 
model this functional form, 

= ^+ 0 ^ 3 / + ^2 J 

The last term models the ‘dip’ in Fig. The location 
of the minimum is given by /o, while a is the overall 
amplitude of the dip and b is the width. We find that 
the frequency location of the dip is very close to the final 
BH’s ringdown frequency, /rd (they agree within our 
uncertainty in calculating /rd), and that the ringdown 
damping frequency /damp is a good approximation to our 
best fit of the width. These quantities are calculated from 
our final mass and spin fits. For these reasons the ansatz 
that we use in practice is. 




(^4:fda 


/damp + (/ - «5 /rd)2 ’ 

(13) 

We find that the parameter 0^5 is in the range 
[0.98,1.04]. The power law terms account for the overall 
trend of the data, and its behaviour at lower frequen¬ 
cies. The constant term translates into a time shift in 
the overall phase, which will be determined by the conti¬ 


nuity requirements of the final IMR phase; see Sec. VIII 
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FIG. 6: Examples of the merger-ringdown (Region Ilb) 

model for three g = 18 configurations where the spin on the 
large BH is xi = {+0-4, 0, —0.8} and three equal-spin g = 1 
configurations (xi ,2 = +0.98,0,-0.95). The configurations 
are ordered top to bottom in the figure. The left panel shows 
the hybrid data, best-fit and final-model predictions over Re¬ 
gion lib. The right panel shows the difference between the 
hybrid data and the best-fit (dashed line) and between the 
hybrid data and the final model (solid line). 


The phase derivative data are fit to Eq. (13) over the 


frequency range [0.45,1.15] /rd- The upper frequency 
1- 15/rd approximates the highest frequency for which 
we have clean NR data. This fitting window was cho¬ 
sen to have some overlap between the intermediate phase 
model, as indicated in Fig. 

The merger-ringdown phase is given by the integral of 
Eq. (pi). 


</>MR =- \ ao + aif - 0:2/ ^ 

V [ 3 


+ 0^4 tan 


-1 f f ~ <^5/rd\ 


/d; 


amp 




(14) 


FIG. 7: The same configurations and layout as in Fig.[^ but 
now showing phase over the intermediate region (Ila). 


2. Region Ila - intermediate 

To bridge the gap between the lowest common fre¬ 
quency of the NR data and the Region lib merger- 
ringdown model, i.e., over the frequency range Mf G 
[0.018,0.5/rd], we use the following ansatz, 

'll <Plnt = A + P2f~^ + • (15) 

The behaviour of the data over this frequency range 
is predominately proportional to 1//. This is not suffi¬ 
cient at higher mass ratios and high anti-aligned spins, 
where /rd can be approximately half that of the equal 
mass non-spinning case. We find that the additional 
f~^ term fits the data well across the entire parameter 
space. The intermediate (Region Ila) ansatz is used over 
the frequency interval [0.018,0.5/rd], but we found that 
the best results were obtained if the data were fit over 
[0.017,0.75/rd]. 

Once again the phase is obtained by integrating 
Eq. 


For the full IMR phase we use the above fit for fre¬ 
quencies larger than 0.5 /rd- At lower frequencies we 
find that 77 0' is fit better by ^ 1/f and we model this 
region (Ila) separately. 

The phase offset that appears as a constant of integra¬ 
tion 0 ^ 0 , and the time-shift term o^i, will both be deter¬ 
mined in the final model by requiring a smooth connec¬ 
tion with the phase from Region Ila. 

Examples of the results are shown in Fig. for six 
configurations at the edges of our calibration parameter 
space. These are equal-spin q = I waveforms with spins 
X = {—0.95,0,0.98} and g = 18 waveforms with spins 
on the larger BH of xi = {—0.8, 0,0.4} (the second BH 
has no spin). In addition to demonstrating that both the 
ansatz and the final model capture the data well, the fig¬ 
ure also illustrates the large differences in the frequency 
range of the merger-rigndown at different points in the 
parameter space. 


4>lnt = I (^/3o + Plf + P2 Log(/) - y / . (16) 

As in Region Hb, the phase-shift due to the constant of 
integration and the time-shift term will be fixed 
by requiring a smooth connection to the Region I phase. 
The results for the corner cases are shown in Fig. 

This completes the modelling of the phase over the 
frequencies for which we have NR data. Region H. We will 
now consider the signal amplitude over the same region, 
before moving on to the inspiral. Region I. 

C. Amplitude 

When we perform the fits to the amplitude across Re¬ 
gion I and Region H, we first factor out the leading order 
PN f~^ behaviour. The resulting data tend to unity 
as the frequency tends to zero, and as with the use of 
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the phase derivative, allows us to identify and model de¬ 
tailed features of the amplitude behaviour; see Fig. 
which shows both amplitude for PN inspiral waveforms, 
and for the full hybrids. 

The normalisation is given by. 


Collocation Point (Mf) 

Value 

Derivative 

fi = 0.014 

Vl — Alns(/l) 

di = A'l^sih) 

/2 = (/l + /3)/2 

V2 = AHyb(/2) 


/s — /peak 

V3 = AyiRifs) 

ds = ^mr(/3) 


lim 

/^o 


y/®^PN(/)' 




2^ 

3 TT^/^ ’ 


and our normalisation factor is therefore. 


TABLE II: Locations of the collocation points, fi, f 2 , f 3 , and 
the corresponding values of the amplitude A(f) and its deriva¬ 
tive A'(f). All information comes from either the inspiral or 
merger-ringodwn models, except for the value V 2 , which is 
read off the input waveform data. 


Aq = 



f-7/e 


(18) 


1. Region Ilb - merger-ring down 


In all previous phenomenological models [8l[l2l[T4], the 
ringdown amplitude has been modelled with a Lorentzian 
function, which is the Fourier transform of the (two- 
sided) exponential decay function. The Fourier trans¬ 
form of the full IMR data instead exhibit an exponential 
decay, as discussed in Paper 1. The amplitude in Region 
lib is fit over the frequency range Mf G [1/1.15,1.2] /rd 
using the following ansatz. 


^MR 

^0 


'Ts/damp 

(/ ~ /rd)^ + (73/damp)^ 


T2(/-/rd) 

^3 f damp 


(19) 


The coefficient 71 G [0.0024, 0.0169] determines the over¬ 
all amplitude of the ringdown. We expect that the fre¬ 
quency width and location of the amplitude peak can be 
inferred from the remnant BH parameters, which moti¬ 
vates the appearance of the ringdown damping frequency 
/damp in Eq. ( p!^ . In practice we find that the width is in¬ 
creased by themctor 73 G [1.25,1.36], and the decay rate 
l/(/damp73) IS modified by the factor 72 G [0.54,1.0339]. 

If we used only the Lorentzian part of Eq. ( p^ , the 
amplitude peak would be located at /rd- With the ad¬ 
ditional exponential factor, the peak is located at 


/peak 


/rd 


/damp73 (^^/^ ^ 


72 


( 20 ) 


2. Region Ila - intermediate 

We now consider the intermediate region (Ila) between 
the end of the inspiral region (I) and the start of the 
merger-ringdown region (lib). 

Fig.| 8 ] shows the TaylorF2 inspiral amplitude in com¬ 
parison to the amplitude in the hybrid data. In some 
cases, we see that we can model the intermediate (Region 
Ila) amplitude by simply smoothly connecting regions I 
and lib. For example, we could fit the four coefficients 
of a third-order polynomial by matching the value of the 
amplitude and its derivative at the end of the Region I 


(nominally Mf = 0.018) and at the beginning of Region 

Eb, /peak- 

In other cases, however, we see that the rescaled ampli¬ 
tude will have a minimum in the intermediate region, and 
a naive connection of the inspiral and merger-ringdown 
regions would not in general locate this minimum cor¬ 
rectly. 

For this reason, we model the intermediate amplitude 
with a fourth-order polynomial. Four of the coefficients 
are fixed (as above), by matching the value and derivative 
of the amplitude at the endpoints of our intermediate 
fit. The lower frequency is chosen as M/i = 0.014, i.e., 
slightly before the end of the inspiral at M f = 0.018, and 
the upper frequency is /s = /peak- The fifth coefficient is 
determined by the value of amplitude of the NR waveform 
at the frequency mid-way between the two, /2 = (/i + 
/3)/2. 

In practice, the amplitude values and derivatives at 
the endpoints are given by the models for Region I and 
Region IIB. The only additional piece of information that 
needs to be modelled from the NR data is the value of 
the amplitude at / 2 . We find that this can be accurately 
modelled across the parameter space by a polynomial 
ansatz in ( 7 , 7 ), as will be described in Sec. |VII[ 

This collocation method is similar to that used in spec¬ 
tral methods. Given an ansatz with n free coefficients we 
require n pieces of information from the data to constrain 
the ansatz and solve the system. In this case we use the 
value of the function at three points, and the derivative 
at two points. The intermediate ansatz is given by 

Eint = Aq ((io + Sif + 62/^ + Ssf^ + ^4/^) 5 ( 21 ) 

and the 6i coefficients are the solution to the system of 
equations. 


Elnt(/l) = ^ 1 , ( 22 ) 

AiM) = V2, (23) 

Alntifs) = ^ 3 , (24) 

4nt(/i) = du (25) 

A'Uh) = ds. (26) 

The frequencies and values are given in Tab. [TTj 

The results of our amplitude model are shown in Figs.[^ 
and which show the same equal-mass and q = IS 
cases as in Fig. The left panels show the full signal 
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Mf 


Mf 


FIG. 8: Hybrid Fourier domain amplitude for three equal mass cases — 0.98, Xi — X 2 = 0 and Xi — X 2 = —0.95, 

indicated by black, orange and green lines respectively. The PN prediction is shown as dashed lines. The left panel shows the 
full fourier domain amplitude, while the right panel shows the fourier domain amplitude but rescaled by Aq^, Eq. (18). 


amplitude, while the right panels show the amplitude 


scaled by the factor, Eq. 18 


The scaled plots indicate that the weakest part of the 
model is that which describes the intermediate Region 
Ila amplitude. This is because the minimum that we 
see in the scaled figures (those in the right panels) is 
captured only through the value of the amplitude at the 
frequency in the middle of Region Ila. If we were in ad¬ 
dition to model the frequency at which the minimum oc¬ 
curs, and prescribe the amplitude value there, the model 
may perform better. We could also, of course, add fur¬ 
ther collocation points. However, we can see from the full 
unsealed amplitude (the left panels) that the amplitude 
is nonetheless very accurately represented, and in addi¬ 
tion, small variations in the amplitude play a far smaller 
role in GW applications (both searches and parameter 
estimation) than the GW phase. 


VI. INSPIRAL MODEL (REGION I) 

We now turn our attention to modelling Region I, i.e., 
the inspiral portion of the waveform, below the frequency 
Mf = 0.018; see Fig. 

The non-spinning m and the first aligned-spin [8] phe¬ 
nomenological models used a PN-like ansatz for the in¬ 
spiral phase, calibrated against PN+NR hybrids. In the 
PhenomC model m, the TaylorF2 phase was used for 
the equivalent of Region I; in that model the inspiral re¬ 
gion ended at O.I/rd- For the parameter space covered 
by our new model, this would corresponds to frequencies 
between Mf ^ 0.005 and Mf ^ 0.012. 

In Paper 1 we presented evidence that the uncalibrated 
SEOBv2 model is currently the inspiral approximant that 
is most consistent with NR data for the inspiral. In 
this section we construct a frequency domain model of 
the SEOBv2 inspiral, up to Mf = 0.018, using our 
SEOBv2+NR hybrids. As discussed previously, we ex¬ 


pect that the SEOBv2 model is sufficiently accurate up 
to this frequency, and very likely to higher frequencies, 
allowing us to match to our merger-ringdown model at 
significantly higher frequencies than was considered rea¬ 
sonable with the TaylorF2 approximant used for Phe¬ 
nomC. 

Note that it is possible, in principle, to cover the pa¬ 
rameter space with an arbitrarily high density of SEOBv2 
waveforms, and use those to calibrate an inspiral model. 
In this paper, however, we use hybrid SEOBv2+NR 
waveforms and therefore calibrate the inspiral model to 
the same points in parameter space as used for the Re¬ 
gion II merger-ringdown models. 


A. Phase 

The inspiral portion Mf G [0.0035,0.018] of the hy¬ 
brids can be accurately modelled with an ansatz con¬ 
sisting of the known TaylorF2 terms for the phase, aug¬ 
mented with the next four higher order PN terms, with 
their coefficients fit to the SEOBv2+NR hybrid data. We 
find that these higher order terms are enough to capture 
the EOB and NR data over this frequency range to a very 
high level of accuracy. 

The full TaylorF2 phase is, 

(f>TF 2 = 27r/ic - (^c - 7r/4 

where (pifB.) are the PN expansion coefficients that are 
functions of the intrinsic binary parameters. Explicit 
expressions are given in Appendix We incorporate 
spin-independent corrections up to 3.5PN order (i = 7) 
[501162], linear spin-orbit corrections up to 3.5PN order 
[63] and quadratic spin corrections up to 2PN order [6^ 
l66] . In re-expanding the PN energy and flux to obtain 
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FIG. 9: Hybrid and model Fourier-domain amplitude for three equal-mass configurations, Xi — X 2 = 0.98, Xi = X 2 = 0 
and Xi = X 2 = —0.95, indicated by black, orange and green lines respectively. The hybrid data are shown by solid lines, and 
the PhenomD model by dashed lines. The left panel shows the full Fourier-domain amplitude, while the right panel shows 
the Fourier-domain amplitude but rescaled by ^ Eq. (18). The short vertical dashed lines mark the three frequency points 
in Tab while the lines at lower and higher frequency coincide with the transition points between regions I and Ila and 
between regions Ila and Ilb respectively. 




FIG. 10: The same quantities as in Fig. but now for three g = 18 configurations, yi = 0.4, X 2 = 0, xi = X 2 = 0 and 
Xi = -0.8, X2 = 0. 


the TaylorF2 phase, we drop all quadratic and higher- 
order spin corrections beyond 2PN order as they would 
constitute incomplete terms in our description. With 
these choices, we are entirely consistent with the current 
state of the LIGO software library [58]. We note that we 
also constructed a full model that incorporated recently 
calculated higher-order terms, specifically quadratic spin 
terms at 3PN order m and cubic spin terms at 3.5PN 
order [68] , but we found no significant difference between 
both constructions. 


Equation (27) includes both spins, xi X 2 , while 
our fit for the coefficients of additional terms will be pa¬ 
rameterized only by x- This means that the final phase 
expression will incorporate some effects from the spins 
of each BH, but, although the model is sufficiently ac¬ 
curate for use in GW astronomy applications across a 
wide range of the two-spin parameter space, it should 
not be considered an accurate representation of two-spin 


effects. We expect the model to be more than sufficient 
for searching for BH binaries with any BH spins within 
the calibration parameter space, or for estimation of the 
parameters (M, 77,x), but we do not recommend its use 
in, for example, theoretical studies of detailed double¬ 
spin effects in binaries. 

The phase ansatz is given by. 


01ns =0TF2(T/'/; S) 

+ “ (^^0 + CTi/ + • 

(28) 

Note that to compute the phenomenological coeffi¬ 
cients the fit is performed over the frequency range 
Mf G [0.0035, 0.019] to achieve an optimal balance be¬ 
tween goodness of fit and accuracy in reproducing phe- 
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FIG. 11: The same analysis as in Figs. |^and|^ but now for 
the inspiral model. 


(q^o, /^O: A) across Region II (See Sec. VIII) are con¬ 

strained analytically; there is only one time and phase- 
shift freedom for the full waveform. This leaves a to¬ 
tal of 17 phenomenological parameters which need to 
be mapped on to the physical parameter space. We 
parametrise the phenomenological coefficients by two 
physical parameters, (77 ,Xpn)- Our model is also depen¬ 
dent on the total mass M of the system through a trivial 
rescaling. 

As in previous phenomenological models nnaini we 
map the phenomenological coefficients in terms of poly¬ 
nomials of the physical parameters, up to second order in 
T] and third order in Xpn, although in this work our poly¬ 
nomial ansatz is expanded around xpn = 1- Note that 
in the fit across the parameter space we use the unsealed 
reduced-spin parameter Xpn, 


nomenological coefficients and to reduce boundary effects 
at the interface between Region I and Region Ila (i.e., 
Mf = 0.018). In practice the fits were performed over 
the (j)' data, as with Region II above. We will see in 
Sec. |IXA| that this model also sufficiently accurately rep¬ 
resents SEOBv2+NR hybrids down to much lower fre¬ 
quencies. 

The results for the three example q = I and q = IS 
configurations are shown in Fig. We see that once 
again our ansatz accurately models the data, and that the 
Fourier-domain phase error is below 0.15 rad for the en¬ 
tire inspiral for the high-mass ratio configurations, while 
for the equal-mass configurations the phase error is typ¬ 
ically an order of magnitude smaller. 

B. Amplitude 


= ^00 + ^ 10 ^ 

(XPN — 1) (Aqi + 

+ (XPN - 1)^ (Ao2 + >^12^ + ^22V^) 

+ (Xpn - 1)^ (Aqs + 5 

where A* indexes both the amplitude and phase coeffi¬ 
cients given by, 

A = { {pj }, {v2 }, {Xj } , }, {Pj }, {oi^j } } . (32) 

V V V V 

S/* V* 

Amplitude Coefficients Phase Coefficients 

The index i G {1,2, 3,4, 5, 6} selects the coefficient vec¬ 
tor for either amplitude or phase and model for Region 
I, Ila or Ilb. Tab. [V| in Appendix [0 contains the values 
of all the mapping coefficients for each phenomenological 
parameter. 


Our model of the inspiral amplitude is based on a re¬ 
expanded PN amplitude, as discussed in Sec. IV of Paper 
1. The base amplitude is given by. 


Apn(/) = Ao'^Aiin/y 


/3 


(29) 


i=0 


where Aq is the leading order / behaviour in Eq. (18). 


The higher order terms that we calibrate are the next 
natural terms in the PN expansion. 


VIII. FULL IMR WAVEFORMS 

By construction, all the regions of the amplitude and 
phase models are joined by C(l)-continuous conditions. 
This ensures the first derivative of the amplitude and 
phase at the boundary between the various regions, which 
are used in analytic calculations, are smooth. We assume 
that this is sufficient and simply join together the piece- 
wise regions with step functions. Our step function is 
defined as 


Alns — A 


PN 




)/3 


(30) 




and, 

VII. MAPPING THE PHENOMENOLOGICAL 
COEFFICIENTS TO PHYSICAL PARAMETERS 


0{f - fo) 


-1, / < fo, 
1, / > fo, 


Xo = ^ [1 ± ^(/ - fo)] ■ 


(33) 


(34) 


Our model has 11 amplitude and 14 phase coeffi¬ 
cients. Howev er four of the amplitude coefficients in Re¬ 
gion Ila (Sec. VC 2) and four of the phase coefficients 


The full IMR phase is determined up to an arbitrary 
time- and phase shift. These shifts are absorbed into 
the constant and linear coefficients of the inspiral part 
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(cro,cri). The constant and linear coefficients of the Re¬ 
gion Ila (ofo, ai) and Ilb models (/3o, Pi) are fixed by the 
requirement of C{1) continuity. 

The full IMR phase is given by the following equation 


^IMr(/) = Plnsif) plntif) 0 0Mr(/) , 




/2 ^ ^/2 ' 


(35) 

where pins is given by Eq. ([28]), pint by Eq. ([^, and puR 
by Eq. (in]), and the transition frequencies are fi = 0.018 
and /2 = 0.5/rd. As noted previously, when evaluating 
the known PN part of pins^ given in Eq. (28), we use the 
full two spin dependence. 

The full IMR amplitude is given by 


Aimr(/) = Ains(/) Aint(/) 0 Amr(/) , 

where Ains is given by Eq. (j^, Aint by Eq. (|2^, and 
Amr by Eq. (19), and where the transition frequencies 
are fi = 0.014 and /2 = /peak, Eq. (20). The amplitude 
is C(l)-continuous by construction. Once again, note 
that the base inspiral PN amplitude includes both spin 
contributions. 

The phase and amplitude coefficients across the (77, x) 
parameter space are shown in Eigs. and We 

see that in general the coefficients vary smoothly across 
the parameter space, and are captured well by our fits. 


IX. MODEL VALIDATION 


To evaluate the accuracy of our model we compute the 
mismatch, defined in Sec. El between the model and 
a set of hybrid waveforms, including the 19 waveforms 
used to calibrate the model (Tab. [Ij), and an additional 
28 waveforms, listed in Tab. EH The additional SpEC 
NR waveforms comprise most of the remaining aligned 
spin simulations in the public SXS catalogue [43]. The 
remaining NR waveforms were produced with BAM. 

In this section we quantify the agreement for each of 
these wave forms against the PhenomD model. We also 
show (Sec. |IXC ) that using additional waveforms in the 
calibration does not significantly change our model, and 
provide evidence that the set of waveforms we have cho¬ 
sen may be close to the minimal set necessary to accu¬ 
rately calibrate our model. 

A further, complementary validation based on time- 
domain transformations is presented in Appendix [Aj 


Mismatches 


In this section we compute the mismatch between Phe¬ 
nomD and all of the hybrid waveforms in Tabs. |T| and III 


The model was calibrated to hybrid waveforms with 
a starting frequency of Mf = 0.0035, but the wave¬ 
forms from many astrophysical compact binaries will be 
detectable by aLIGO and AdV from much lower frequen¬ 
cies. We assume that the minimum mass for one of the 


# 

Code/ID 

Q 

Xi 

X2 

B1 

SXS:BBH:0159 

1 . 

-0.9 

-0.9 

B2 

SXS:BBH:0154 

1 . 

-0.8 

-0.8 

B3 

SXS:BBH:0148 

1 . 

-0.438 

-0.438 

B4 

SXS:BBH:0149 

1 . 

-0.2 

-0.2 

B5 

SXS:BBH:0150 

1 . 

0.2 

0.2 

B6 

SXS:BBH:0170 

1 . 

0.437 

0.437 

B7 

SXS:BBH:0155 

1 . 

0.8 

0.8 

B8 

SXS:BBH:0153 

1 . 

0.85 

0.85 

B9 

SXS:BBH:0160 

1 . 

0.9 

0.9 

BIO 

SXS:BBH:0157 

1 . 

0.95 

0.95 

Bll 

SXS:BBH:0158 

1 . 

0.97 

0.97 

B12 

SXS:BBH:0014 

1.5 

-0.5 

0 . 

B13 

SXS:BBH:0008 

1.5 

0 . 

0 . 

B14 

SXS:BBH:0013 

1.5 

0.5 

0 . 

B15 

SXS:BBH:0169 

2 . 

0 . 

0 . 

B16 

BAM 

2 . 

0.5 

0.5 

B17 

BAM 

2 . 

0.75 

0.75 

B18 

BAM 

3. 

-0.5 

-0.5 

B19 

SXS:BBH:0036 

3. 

-0.5 

0 . 

B20 

SXS:BBH:0168 

3. 

0 . 

0 . 

B21 

SXS:BBH:0045 

3. 

0.5 

-0.5 

B22 

SXS:BBH:0031 

3. 

0.5 

0 . 

B23 

SXS:BBH:0047 

3. 

0.5 

0.5 

B24 

BAM 

4. 

-0.25 

-0.25 

B25 

BAM 

4. 

0.25 

0.25 

B26 

SXS:BBH:0060 

5. 

-0.5 

0 . 

B27 

SXS:BBH:0056 

5. 

0 . 

0 . 

B28 

SXS:BBH:0166 

6 . 

0 . 

0 . 

B29 

BAM 

10 . 

0 . 

0 . 


TABLE III: Additional waveforms used to verify the model, 
but not used in its calibration. 


compact objects is given by the typical mass of a neu¬ 
tron star i.e., Mns ^ 1.4M0. The total mass of the 
binary can then be no lower than Mmin = (^ + l)ArNS 
for configurations with mass-ratio q. Our goal is to pro¬ 
duce a model that is accurate for binaries that can be 
detected from 10 Hz down to either 12 Mq ISO], or Mmin, 
if this exceeds 12 Mq, which is the case for systems with 
q ^ S. At 10 Hz, the waveform frequency of a 12 Mq bi¬ 
nary is Mf ^ 0.0006, and so in this section we compare 
our model to much longer hybrids that extend down to 
Mf = 0.0006. 

The results are presented in Eig. The left panel 
uses the aLIGO design sensitivity zero-detuned, high- 
laser-power noise curve with [/min,/max] = [10, 8000] Hz 
|69| . The worst mismatch is for the {q, Xi, X2} = {6, 0, 0} 
at low masses which tends towards a mismatch of 3% at 
12 Mq. All other mismatches fall below 1% with the ma¬ 
jority distributed around 0.1%. We note, however, that 
the fitting factors for the waveforms in the model (i.e., 
matches optimised over binary parameters, which is the 
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FIG. 12: Phase coefficients for region I and II. The calibration points and the model, extrapolated to the boundary of the 
physical parameter space are shown. 



FIG. 13: Amplitude coefficients for region I and Ilb. The calibration points and the model, extrapolated to the boundary of 
the physical parameter space are shown. 
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FIG. 14: Intermediate (Region Ilb) amplitude coefficient. 
The calibration points and the model, extrapolated to the 
boundary of the physical parameter space are shown. 


relevant quantity for searches) are better than 0.999 in all 
cases we have considered. In particular, at low masses the 
mismatch between different options of inspiral approxi- 
mant will be much larger than the mismatch between 
PhenomD and our hybrid waveforms; the dominant error 
is in our uncertainty of the true inspiral waveform, and 
not in our model; this will be made clearer in Sec. 

The right panel in Fig 15 shows the same calculation 
but using the predicted noise curve for early aLIGO sci¬ 
ence runs [41], with a lower frequency cut-off of 30 Hz. 
Due to the change in shape of the noise curve and lower 
frequency cut-off the mismatches improve such that all 
mismatches are comfortably below 1%. This gives a more 
realistic idea of the performance of our model during the 
initial science run of the advanced detectors. 

In both panels, the highlighted cases are those at the 
edges of the calibration region of parameter space. We 
note that the worst mismatches are for high mass ratios 
and large spins. This suggests the region of parameter 
space that requires the most improvement in future mod¬ 
els — although it is clear that for all of these configura¬ 
tions the model is well within the accuracy requirements 
for the second-generation detectors. 


B. The effective spin approximation 

The phenomenological fits to the waveform phase and 
amplitude are parameterised by the weighted reduced 
spin, X, Eq. Q. This is an approximation, based on the 
observation that the dominant spin effect on the inspi¬ 
ral phase is due to this combination of the two spins, xi 
and X 2 - This approximation is not expected to be valid 
through the merger and ringdown; in the ringdown the 
waveforms will be characterized by the final spin. The 
model was produced using mostly equal-spin xi = X 2 
waveforms, and in general may not be accurate for sys¬ 
tems with unequal spins. 

However, we have seen in the previous Sec. |IX A| that 
our model agrees well with all available hybrid wave¬ 
forms, including several with unequal spins. This in¬ 


cluded only four unequal-spin configurations that were 
not included in the calibration, and none were high- 
aligned-spin systems. 

We expect that the reduced-spin approximation will 
perform worst for high mass ratios and high aligned spins. 
If we consider pure PN inspiral waveforms, we find, for 
example, that a system with mass-ratio 1:3 and total 
mass of 12 Mq, with = 1 and X 2 = —1, that the 
match against the corresponding reduced-spin waveform 
(with X = 0.655) is less than 0.8. However, if we consider 
a configuration where the larger BH has an anti-aligned 
spin, xi = ~lW 2 = 1, then the match with the corre¬ 
sponding reduced-spin waveform (x = —0.655) is much 
better, 0.955. 

This example was only an illustration. The perfor¬ 
mance of the reduced-spin approximation at low masses 
does not concern us in the PhenomD model, where we 
use both spins xi and X 2 to generate the base TaylorF2 
phase. What we wish to know is how well the approxima¬ 
tion holds for high-mass systems, where the late inspiral, 
merger and ringdown dominate the SNR. Those systems 
are described by our merger-ringdown Region H model, 
for which the spin dependence is parameterized only with 
X- 

We have produced one high mass-ratio, high-spin NR 
simulation to compare with, q = S and xi = X 2 = 0. 
Fig. shows the mismatch between this hybrid wave¬ 
form and the PhenomD model. As we expect in this 
region of the parameter space, the poor quality of the 
reduced-spin approximation causes a mismatch that ex¬ 
ceeds our 1% threshold for all masses. However, if we 
calculate fitting factors (i.e., minimise the mismatch with 
respect to the model parameters (r^, x), as done in a GW 
search and, effectively, in parameter estimation), then we 
find deviations from unity of below 0.05% for all masses. 
We also find biasses of less than 1% in the total mass, 
less than 2% in the symmetric mass ratio, and less than 
0.005 in the reduced spin, x- We expect these biasses to 
be far less than the statistical uncertainties in these quan¬ 
tities for observations with second-generation detectors, 
and so we conclude that the reduced-spin approximation 
will not impose any limit on the science potential of these 
detectors. 

Studies with the SEOBNRv2 model support this con¬ 
clusion. Although we do not expect that model to be 
accurate through the merger-ringdown for high spins, as 
we will see in Sec. it is likely that its qualitative be¬ 
haviour with respect to parameter variations is approxi¬ 
mately correct, and the model allows us to study the be¬ 
haviour of the reduced-spin approximation over the entire 
calibration parameter space of our model. 

Although the reduced-spin approximation will not 
limit our ability to measure x^ could argue that 
it nonetheless prevents any measurement of individual 
spins. We have argued in previous work [34| that it may 
be difficult to measure both BH spins even if we have 
a double-spin model. A further study, which provides 
much stronger evidence for this claim, will be published 
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FIG. 15: Mismatches of the PhenomD model against all 48 available hybrid waveforms. The highlighted configurations are 
those closest to the edge of the (?7,x) parameter space as well as the case with the worst mismatch (g,Xi,X 2 ) = (6,0.,0). The 
majority of cases show mismatches well below 1%. Left: Mismatches using the aLIGO design sensitivity noise curve (zdethp) 
with a lower frequency cut off of 10 Hz. Right: Early aLIGO noise curve with a 30 Hz cut off. 



FIG. 16: Mismatch between a g = 8 , xi = 0-8, X 2 = 0 
SEOBv2+NR hybrid, and the PhenomD model. We see that 
the mismatch exceeds our 1% threshold everywhere. How¬ 
ever, the fitting factor is everywhere better than 0.9995, with 
negligible parameter biases (see text). 


in the near future [35]. In practice the measurable intrin¬ 
sic parameters of the binary will be (M, 77,x), and these 
are the parameters of our model. 


C. Calibration Set of waveforms 

The construction of previous phenomenological mod¬ 
els imiHii] suggested that the parameter dependence 
of the coefficients in our models depend sufficiently 
smoothly across the parameter space that each coeffi¬ 
cient can be presented by a low-order polynomial in each 
parameter, and therefore we require only 4-5 waveforms 
for each direction in parameter space. This expectation 
is borne out in the current model, where we use four val¬ 
ues of the mass ratio (1, 4, 8 and 18) and four or five 
values of the spin at each mass ratio. 

In this section we consider versions of the model con¬ 
structed with more (or less) calibration waveforms. We 
find that our small set of 19 calibration waveforms is 
just as accurate as a model that is calibrated against a 
much larger set of 48 waveforms. To quantify this test we 


compute the maximum mismatch of four distinct models 
against all hybrid waveforms used in this paper, i.e., the 
48 waveforms in Tabs. [Tj and m 

Fig indicates four choices of parameter-space cov¬ 
erage. The first set is the largest, and includes all 48 
configurations indicated in the figure. The second set in¬ 
cludes 25 waveforms, but only at mass ratios 1, 4 8 and 
18, and does not include all available spin values at mass 
ratios 1 and 8. The third set consists of the 19 waveforms 
that we use for our final model. The fourth set is more 
sparsely sampled in spins, with only three spin values at 
each mass ratio, and only 12 waveforms in total. 


Four models were constructed, each using the same 
prescription, except for the Set-4 model, for which we 
used a lower-order fit in the x direction, since in general 
we cannot expect to fit four coefficients with only three 
spin values. 


The results are summarized in Tab. nvi We calculate 
the mismatch between each of the models and all 48 hy¬ 
brids, over the same mass range used in Sec. |IX A] using 
the early aLIGO noise curve with a 30 Hz cut off. For 
each hybrid we calculate the largest mismatch in that 
mass range. The table indicates the number of configu¬ 
rations for which we find mismatches larger than 0.1%, 
1% and 3% for each model. As we have already seen 
in Fig. the fiducial Set-3 model has mismatches of 
less than 1% for all configurations. We find that in¬ 
creasing the number of calibration waveforms does not 
significantly improve the model’s performance. 


We also see that if we further reduce the number of 
calibration waveforms, as in the Set-4 model, then the 
accuracy of the model drops significantly. For this model 
there are now three configuration with mismatches worse 
than 1%, and one configuration with a mismatch worse 
than 7%. We therefore conclude that, in the sense of 
the simple comparison that has been performed here, the 
Set-3 model represents the optimal choice of calibration 
waveforms. 
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18 8 4 1 



n 

FIG. 17: Four sets of calibration waveforms. Set 1 (48 wave¬ 
forms) is indicated in red, Set 2 (25 waveforms) in green, Set 3 
(19 waveforms, used for the final PhenomD model) in orange, 
and Set 4 (12 waveforms) in blue. 


Model 

# waveforms >0.1% 

> 1% 

> 3% 

max Ad (%) 

Set 1 

48 

19 

0 

0 

0.94 

Set 2 

25 

27 

0 

0 

0.83 

Set 3 ( 

*) 19 

29 

0 

0 

0.87 

Set 4 

12 

37 

3 

1 

7.82 


TABLE IV: Comparison of models constructed with differ¬ 
ent sets of calibration waveforms. The table shows, for each 
calibration set (see Fig. 17), the number of verifications wave¬ 
forms (out of 48) for which there is a mismatch M above 0.1%, 
1%, or 3%, over the same mass range used in Sec. |IX Ab using 
the early aLIGO noise curve with a 30 Hz cut off. We see 
that with a small set of 19 waveforms we achieve comparable 
mismatches to models which used larger sets of calibration 
waveforms, and that using less waveforms significantly de¬ 
grades the quality of the model. Set 3 is used for the final 
model. 


X. MODEL VS MODEL COMPARISONS 

We have demonstrated the high degree of fidelity of 
PhenomD to both the waveforms that were used in cali¬ 
brating the model and to those that were not. Without 
further comparisons to NR waveforms we cannot rigor¬ 
ously quantify the accuracy of our, or indeed any, wave¬ 
form model. However, it is reasonable to assume that if 
two independent waveform models agree over a portion of 
the parameter space then we can gain some well-founded 
confidence in their accuracy. 

The computational cost of the SEOBNRv2 model 
makes it difficult to make detailed comparisons across 
the entire parameter space with high resolution in (r^, x). 
However, based upon the recent work in Ref. m, 
a reduced order model (ROM) of SEOBNRv2, called 
SEOBNRv2_ROM, has been developed [71]. This is a 
fast, frequency-domain approximation to the SEOBNRv2 
model that has a worst mismatch against SEOBNRv2 of 
1%, but in general mismatches are better than ^ 0.1%. 


SEOBNRv2_ROM is a two spin model which can be 
used to estimate SEOBNRv2 waveforms with symmet¬ 
ric mass-ratios r] G [0.01,0.25] and spins Xi ^ [—1,0.99]. 
The ROM can be used over the frequency range Mf G 
[0.0001, 0.3] Note that the underlying SEOBNRv2 model 
was calibrated to NR waveforms up to mass-ratios 1:8 
and spins up to 0.5 (except along the equal mass line 
where spins in the range [—0.95,0.98] were used). The 
merger-ringdown parts of the PhenomD and SEOBNRv2 
models are almost completely independent of one another 
with the only common features being that they share 
some of the same calibration waveforms, i.e., the ones 
from the public SXS catalogue and also the same under¬ 
lying EOB Hamiltonian. 

During the following comparison we restrict the com¬ 
putation of the mismatch to the frequencies of the SEOB- 
NRv2_ROM, namely [0.0006,0.135], using the design 
sensitivity noise curve with a lower frequency cut-off of 
10 Hz as in previous sections. 

We noted earlier that the PhenomD model is modu¬ 
lar, and we can use alternative models of either the in¬ 
spiral or merger-ringdown regions as we wish. In the 
following comparisons we consider three versions of the 
model. One is the full PhenomD model that we have 
presented in the previous sections. In comparisons with 
SEOBNRv2_ROM at low masses, the mismatch is dom¬ 
inated by differences between the uncalibrated SEOBv2 
model that we used to calibrate the inspiral of PhenomD, 
and the calibrated SEOBNRv2 model; it is a reflection 
of a different choice of inspiral approximant, and not 
the inherent accuracy of either model. For this reason 
we also perform a second set of comparisons, where we 
use SEOBNRv2_ROM for the inspiral (Region I) part 
of PhenomD; the merger-ringdown (Region II) remains 
unchanged. This allows us to compare PhenomD and 
SEOBNRv2_ROM over only the merger-ringdown, and 
also illustrates the flexibility of the PhenomD model in 
using alternative inspiral approximants. Finally, we re¬ 
place the inspiral part of the PhenomD with TaylorF2. 

The results of our comparions are shown in Fig. 
Each panel shows the mismatch in percentage between 
the PhenomD model and SEOBNRv2_ROM (left 
column) and between [SEOBNRv2_ROM-inspiral + 
PhenomD-merger-ringdown] and SEOBNRv2_ROM 
(middle column), and between [TaylorF2-inspiral + 
PhenomD-merger-ringdown] and SEOBNRv2_ROM 
(right column). The calculations were performed over 
mass ratios [1,100], spins in the range [—1,0.99] and 
for the total masses [12, 20, 50,100,150] Mq. Overlaid in 
white dots are the calibration points of the PhenomD 
model. It is instructive when studying these plots 
to recall that the common region of parameter space 
calibration is up to mass-ratios 1:8 {r] ^ 0.01) and spin 
[—0.5, 0.5], except along the equal mass line where the 
spins range from [—0.95,0.98]. 

We focus first on the low-mass configurations (M < 
50Mq). We see that the agreement between PhenomD 
and SEOBNRv2_ROM is in general quite poor — some 
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FIG. 18: Mismatch comparisons between the SEOBNRv2_ROM model, and three versions of PhenomD. Left: the final 
PhenomD model. Middle: SEOBNRv2_ROM is used for the inspiral part of PhenomD, i.e., up to Mf = 0.018. Right: 
TaylorE2 is used for the inspiral part of PhenomD. See text for discussion. 


parts of the common calibration region of both models 
show mismatches greater than 3%, e.g, for anti-aligned 
spins. This is not necessarily due to the innaccuracy of 
either model. We have seen in Fig. that PhenomD 
typically has matches of better than 1% against our hy¬ 
brid waveforms, which demonstrates that the model ac¬ 
curately reproduces the uncalihrated SEOBv2 model at 
low frequencies. Therefore, we expect that the poor mis¬ 
matches between PhenomD and SEOBNRv2_ROM at 
low masses are due to differences between SEOBv2 and 


the calibrated SEOBNRv2 inspiral. This expectation is 
borne out in the middle column, where the SEOBv2- 
based PhenomD inspiral is replaced with the SEOB- 
NRv2_ROM inspiral. Now the modified PhenomD and 
SEOBNRv2_ROM models differ only in their descrip¬ 
tion of the merger-ringdown, and should agree well at 
very low masses, where the merger-ringdown contributes 
little signal-to-noise ratio (SNR). This is what we find: 
at 12 Mq the mismatches are better than 1% for most 
of the parameter space. The merger-ringdown still has 
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some influence, increasing the mismatches for high-spin 
and high-mass-ratio systems, but in general the agree¬ 
ment is extremely good. 

Although the uncalibrated SEOBv2 and the calibrated 
SEOBNRv2 inspirals show poor matches at low masses, 
we note that both are still consistent with our full NR 
data at higher frequencies, and both are adequate op¬ 
tions for an inspiral description, as we discussed in de¬ 
tail in Paper 1, and also in Sec. |IV| above. The right 
panel illustrates how the model would change if we in¬ 
stead used TaylorE2 for the inspiral. At the matching fre¬ 
quency with the merger-ringdown model {Mf = 0.018) 
the TaylorE2 phase disagree (in the sense of the time- 
shift analysis in Paper 1) at a level that makes it diffi¬ 
cult to smoothly connect them over large regions of the 
parameter space. This, in addition to the differences be¬ 
tween TaylorE2 and SEOBv2(NR) at low frequencies, in¬ 
troduces high overlaps over all but a small strip of pa¬ 
rameter space. 

As we progress down the table of plots to higher 
masses, the merger-ringdown contributes more power to 
the SNR, and the results of the left and middle compar¬ 
isons agree more. At 150 Mq, where the contribution 
from the inspiral (taken here as Mf < 0.018) is negligi¬ 
ble, we see that the two comparisons are almost identical. 
The poor agreement between TaylorE2 and our merger- 
ringdown model at Mf = 0.018 continues to lead to large 
mismatches. 

We now focus on the high-mass configurations (M > 
50Mq), and the left panels that directly compare Phe- 
nomD and SEOBNRv2_ROM. It is evident that the re¬ 
gion of agreement between the two models follows closely 
the region of common calibration points. Indeed, it is 
very encouraging that there is a high level of agreement 
between these two independent models even up to high 
mass-ratios of 1:18 and towards large negative spin val¬ 
ues. 

The positive spin section shows a different behaviour. 
At high masses (i.e., where the merger and ringdown are 
in the detector’s most sensitive frequency range), there is 
a sudden drop in the agreement between the two models 
at mass-ratios larger than equal mass and spin greater 
than ^ 0.75. 

PhenomD is calibrated to two high-spin unequal- 
mass cases, (g, Xi,X 2 ) = {(4,0.75,0.75), (8,0.85,0.85)}, 
and we have one additional case for verification, 
(2,0.75,0.75). These are the waveforms AlO and A15 
from Table [T] and B17 from Table m respectively. As 
we have already seen, PhenomD has better than 1% mis¬ 
match to all theses cases and therefore the poor mis¬ 
matches are unlikely due to errors in PhenomD. We note 
that these cases are well outside the calibration region 
of the SEOBNRv2 model, and we therefore suspect that 
the accuracy of its description of the merger-ringdown 
degrades significantly for high spins. 

Our results also suggest that, despite the lack of cali¬ 
bration waveforms at high anti-aligned spins, the SEOB- 
NRv2 model remains accurate in that region of parameter 



FIG. 19: Mismatch of PhenomD (solid) or SEOB- 

NRv2_ROM (dashed) against cases AlO (orange), A15 
(green) and B17 (black). 


space, and the relatively good agreement between the two 
models even for nearly extreme anti-aligned spins sug¬ 
gests that additional calibration waveforms, while they 
would be valuable, are less crucial in those cases. 

We also observe poor mismatches for very high mass 
ratios. However, since this is outside the calibration re¬ 
gion of both models, we cannot conclude which (if either) 
is correct. 

To illustrate further the disagreement between Phe¬ 
nomD and SEOBNRv2 at unequal masses and high spins, 
we consider in more detail the three NR configurations 
that we have available. Eig. shows mismatches be¬ 
tween pure NR waveforms {not the hybrids) for each of 
these cases, and against the PhenomD and SEOBNRv2 
models, using the techniques discussed in Sec. El The 
mismatch against SEOBNRv2 is above 1% for all masses, 
and can be as high as 10%. We have reproduced these 
plots using SEOBNRv2 waveforms generated from the 
LAL code, and the results are indistinguishable; the poor 
mismatches cannot be attributed to any errors in the 
ROM construction. 

We therefore conclude that the merger and ringdown 
are not accurately represented in the SEOBNRv2 model 
for high spins. This does not detract from the power of 
the EOB NR approach, but simply illustrates that we 
should not expect any merger-ringdown model to be ac¬ 
curate outside its region of NR calibration. The same 
applies to our PhenomD model; we cannot make any 
statements on its accuracy for spins with x ^ 0.85, other 
than for equal-mass systems. 


XI. SUMMARY AND DISCUSSION 

We have presented a new phenomenological model of 
the GW signal from the inspiral, merger and ringdown 
of aligned-spin BH binaries, PhenomD. The new model 
is calibrated to hybrid EOB+NR waveforms that cover 
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the largest region of parameter space of any aligned-spin 
model to date — mass ratios up to 1:18 and spins up 
to ajm ^ 0.85. The inspiral and merger-ringdown are 
described by three separate models, allowing high accu¬ 
racy over the full frequency range detectable by aLIGO 
and AdV, and also making the model modular: the in¬ 
spiral and merger-ringdown parts can easily be modified 
or replaced if improved or extended models (e.g., to a yet 
larger region of parameter space) become available. 

The inspiral part of our hybrids consists of uncalihrated 
SEOBv2 waveforms. We have shown in Paper 1 that 
the SEOBv2 waveforms are the most consistent with our 
NR simulations over the full parameter space that we 
consider, and we choose to use uncalibrated SEOBv2 to 
produce a model that is fully independent of the NR cal¬ 
ibration done to produce SEOBNRv2. 

The merger-ringdown part of the hybrids (i.e., the NR 
waveforms) have a common lowest frequency of Mf ^ 
0.018, and so this is the frequency at which we switch 
from the inspiral to the merger-ringdown model. 

The final model has mismatches against both the 19 
calibration hybrids and an additional 29 verification hy¬ 
brids, of typically better than 1% for all masses. The 
mismatches are shown in Eig. and demonstrate that 
we have faithfully modelled this region of the aligned-spin 
parameter space. 

The model is parameterized by the binary’s symmet¬ 
ric mass ratio, r^, and a normalized reduced effective spin 
parameter, y, defined in Eq. Q- A parameterization 
in terms of a weighted sum of the two BH spins has 
been used in previous Phenom models EHH, and is mo¬ 
tivated by the leading-order spin effect on the inspiral 
phasing [33 El EH, and demonstrations of its efficacy 
for merger-ringdown [34]. In this paper we show that the 
reduced-spin approximation becomes inaccurate only for 
high-spin unequal-mass systems, but in these configura¬ 
tions the parameter errors due to our approximation ap¬ 
pear to be smaller than statistical errors in the spin and 
mass-ratio measurements with aLIGO and AdV. This 
implies that it will be difficult to measure both BH spins 
in GW measurements; this will be considered in more 
detail in a forthcoming paper |35| . 

We have compared the new PhenomD model with the 
state-of-the-art SEOBNRv2 model, and found that the 
two models agree well over their common region of cali¬ 
bration, which is mass ratios up to 1:8, and spins up to 
a/m ^ 0.5 (and near-extremal spins for equal-mass sys¬ 
tems). At low masses the agreement is not good, but we 
show that this is due to differences between the calibrated 
and uncalibrated SEOBv2 inspiral descriptions. 

Outside the common calibration region, the two mod¬ 
els show significant disagreement, in terms of their mis¬ 
match. This is particularly true for high aligned spins. 
Given that PhenomD was calibrated to several high-spin 
unequal-mass simulations (spins of 0.75 or 0.85), while 
SEOBNRv2 was calibrated to spins of no higher than 
0.5 for unequal-mass configurations, we conclude that 
SEOBNRv2 does not accurately capture the merger and 


ringdown for these systems. We expect, however, that its 
performance will become comparable to PhenomD when 
calibrated to additional NR waveforms. 

The broader conclusion we draw from these results is 
that high-aligned-spin systems deserve greater attention 
in future modelling efforts. The PhenomD model was 
calibrated to only two high-aligned-spin binaries, but it 
is clear that a larger number of NR simulations in this 
region of parameter space will benefit GW astronomy. 

The PhenomD model involves 17 coefficients that are 
mapped across the parameter space with polynomials up 
to second order in r] and up to third order in y. Al¬ 
though the total number of coefficients is similar to the 
previous PhenomC model, the development of a refined 
ansatz for each frequency region allows us to more accu¬ 
rately model a wider range of features of the waveforms. 
This is described in more detail in Paper 1. We have also 
carefully tuned each ansatz, and our parameter-space fits, 
to ensure that the model produces physically reasonable 
results outside the calibration region, and that the wave¬ 
forms show no pathological features when converted to 
the time-domain (Appendix]^. These modifications sig¬ 
nificantly improve the model beyond previous Phemom 
models, in addition to increasing the range of calibration 
and lowering the mismatch error. 

In previous work we have shown that models for 
generic (precessing) binaries can be produced by “twist¬ 
ing up” an aligned-spin model. The PhenomP model 
exploits that idea, but to date has been based on the 
PhenomC model, which limits its applicability to mass 
ratios ^ ^ 4. With the advent of PhenomD, we will be 
able to make PhenomP valid to much higher mass ratios 
and higher values of the parallel component of the spin. 
This simple replacement of PhenomC with PhenomD in 
the LIGO-Virgo LAL code has already been tested, and 
will be made available in the near future. 
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Appendix A: Time-domain conversion 


Our PhenomD model is formulated entirely in the fre¬ 
quency domain, which is a great advantage for perform¬ 
ing fast GW searches and parameter estimation studies. 
However, our construction process started with data in 
the time domain, and physical signals are smooth func¬ 
tions in both the frequency and time domain. Therefore, 
it is desirable to check how our model transforms from 
the frequency domain back into the time domain via a 
straightforward inverse Fourier transformation. 

This serves also as an independent, powerful sanity 
check. The previous PhenomC model m, for instance, 
quickly develops a pathological behavior in the time do¬ 
main once the parameters leave the calibration region, 
which is a result of steep transitions caused by extrapo¬ 
lating fitting coefficients. We do not find these features 
for our new PhenomD model. 

Before applying the inverse Fourier transformation, we 
multiply our model with a variant of the Planck taper 
function 


TU) 


'o, 

< [exp(^ + =^)+l 




-1 


/</l 

/l < / < /2 > 
/ > /2 


where /2 is the smallest frequency that we want to repre¬ 
sent in the time-domain data (which become infinitively 
long for /2 ^0). In order to avoid a sharp transition, 
which would introduce unphysical oscillations, T uses an 
extra cushion, / G (/i,/2), in which the frequency do¬ 
main amplitude smoothly increases from zero to their 
correct value. We typically set fi = O.8/2. 

We perform the Fourier transformation numerically, 
which requires us to define a suitable sampling rate in 
the time and frequency domain. From our model, we 
find that the amplitude has dropped several orders of 
magnitude for frequencies Mf > 0.25, so we can choose 
any sampling with At/M < 2 which in turn is solely 
determined by the largest frequency we include in our 
frequency-domain data. 

The frequency-domain sampling, on the other hand, is 
determined by the total length of the signal in the time 
domain, which is information we do not have a priori ac¬ 
cess to. However, in the spirit of the stationary-phase ap¬ 
proximation that typically relates the time-domain phase 
derivative to the frequency {d(j)(t)/dt ~ 27r/), we approx¬ 
imate 
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In (A2), we have introduced an extra factor of 1/2 to 
account for the negative-frequency content of real-valued 
signals (just like in the usual sampling theorem), and 


when choosing A/ we usually apply another factor of 
1/2 as safety margin. 

The time-domain waveforms we obtain this way can be 
compared to the original NR data, and for corners of the 
parameter space used for calibration we show the results 
in Fig.[^ Note that a small overall time and phase shift 
was applied to the model, as these parameters are not 
meant to faithfully capture the arbitrary choices made in 
the original NR simulations. No other optimization has 
been applied. The agreement visible in Fig. [^through¬ 
out the late inspiral, merger and ringdown is remarkable 
and a strong i ndicati on (in additional to the matches pre¬ 
sented in Sec. IX Aj ) that our hybridization, fitting and 
interpolation procedures accurately represent the original 
data. 

In addition to complementing the model validation, we 
may also use the time-domain representations as a visual 
sanity check, even outside the model’s calibration region. 
As mentioned above, this proved to be a powerful test 
of the previous PhenomC model that failed to produce 
reasonable time-domain waveforms in many parts of the 
parameter space outside its calibration range. PhenomD, 
however, does not show any pathological behavior outside 
its calibration region, neither in the time nor frequency 
domain. We illustrate this fact in Fig. [^by showing a 
case where the model parameters have been extrapolated 
to mass ratio 50 and near-extremal spins Xi = X2 = 0.99. 
While such a plot is by no means a guarantee that the 
waveforms are accurate in this regions of the parame¬ 
ter space, it is reassuring that our new model is much 
more robust in its extrapolation, which will allow GW 
search algorithms to use our model slightly outside its 
calibration region, even if we cannot vouch for the level 
of accuracy there. 


Appendix B: PN coefficients 


For the convenience of the reader, we list below the PN 
coefficients implemented in our model. We incorporated 
spin-independent corrections up to 3.5PN order (i = 7) 
[501162], linear spin-orbit corrections up to 3.5PN order 
[63] and quadratic spin corrections up to 2PN order [64]- 
l66] . Our re-expansion strategy follows the choices made 
in the current state of the LIGO software library [58] as 
discussed in Sec. ErAi 


Following (27), we express the frequency-domain phase 


as 


^tf 2 = 27r/4 - (/9c - 7r/4 

The individual masses and spin parameters, rui and Xi 
(i = 1,2), are encoded in the following parameter combi¬ 
nations. 


M = mi + m2 
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FIG. 20: Time-domain PhenomD waveforms (solid, light blue online) and NR waveforms (dashed, red online) for corners of the 
parameter space used for calibration. We plot the plus polarization normalized by the extraction radius, and the binary’s 
parameters are indicated by the mass ratio q = milm 2 and the two spin parameters Xi,X 2 - 
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FIG. 21: Time-domain representation of the PhenomD model outside its calibration region, here for mass ratio 50 and spin 
parameters of xi — X 2 = 0.99. 


r] = mim2/M^, (B2) 

6 = (mi — m2)/M, (B3) 

Xs = (xi + X2)/2, (B4) 

Xa = (Xi - X2)/2. (B5) 

The expansion coefficients are then given by 

ipo = 1, (B6) 


<^i = 0, 


<^2 


3715 


ips = -IOtt + 


557/ 

113^Xa 

3 


113 767/\ 

“3 ^ J 


(B7) 

(B8) 

(B9) 


ip4, = 


15293365 27145r? 30857?2 


508032 

If 5 = [H-log(7rM/)] 

_ 11583231236531 
“ 4694215680 

6848 


504 72 

386457r Q5 tt7 ] 


405 


405 


5“ + ^<^XaXs + ( 5“ + TT ) Xsi 


756 9 

68487^; 6407r^ 


6 - 


732985 


1407/\ 


21 


2268 9 

15737765635 


3048192 


Xa - 

22557r2 

12 


405 57/ 

732985 242607/ 3407/^ 

2268 ^ 81 9 

760557/2 1278257/3 


7/- 


1728 


1296 


2270 

log(647rM/) H- T^TrSxa + 


/22707r 


(P7 — 


63 ov .y 3 

770966757r 3785157r7/ 740457r7/2 


V 3 


5207r7/ Xs, 


254016 


1512 


756 


6 - 


25150083775 268049357/ 19857/^ 


3048192 


6048 


48 


Xa 


25150083775 105666555957/ 10421657/2 53457/' 


3048192 


762048 


3024 


36 


Xs- 


Xs 


(BIO) 

(Bll) 


(B12) 


(B13) 


As discussed in Sec. |VIB] and Sec. IV in Paper 1, our 
inspiral amplitude model is based on a re-expanded PN 
amplitude. The expansion coefficients of Eq. (29) are 
given by 
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Appendix C: Phenomenological Coefficients 


calculated under the parametrization ( 77 , Xpn)- 


The values of the coefficients for the mapping functions 
given in Eq. (31) are shown in Tab. [V[ These values are 
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